function [effect,effect_p,spike_samples,slab_samples,llik_samples] = sample_spike_slab_logit_reg(Y,X,prior_odds,tausq,burnin,ngibbs)
%
% [effect,effect_p,spike,slab,llik] = sample_spike_slab_logit_reg(y,X,prior_odds,tausq,burnin,ngibbs)
%
%  y          = n x r response vector (0, 1)
%  X          = n x p predictor matrix
%  prior_odds = r x p prior log odds ratio
%  tausq      = 1 x 1 prior variance
%
%  effect     = r x p
%  effect_p   = r x p
%  spike      = r x p x ngibbs
%  slab       = r x p x ngibbs
%  llik       = 1 x (ngibbs + burnin)
%
% code: Yongjin Park, ypp@csail.mit.edu
%

    [n,~] = size(X);
    [~,r] = size(Y);

    % add intercept
    X = [X, ones(n,1)];

    prior_odds = [prior_odds, zeros(r,1)];

    [~,p] = size(X);

    assert(numel(prior_odds) == r*p,'prior odds for each spike element');

    nslice = 1;
    slice_width = 1;

    spike = ones(r,p,'single');
    slab = zeros(r,p,'single');

    spike_samples = ones(r,p,ngibbs,'single');
    slab_samples = zeros(r,p,ngibbs,'single');
    llik_samples = zeros(1,ngibbs+burnin,'single');

    for iter = 1:(burnin + ngibbs)

        % sample slab by slice sampling
        log_dist_slab = @(slab) sum(sum(Y.*(X*bsxfun(@times, slab, spike)') - log(1+exp(X*bsxfun(@times, slab, spike)')))) ...
            - 0.5*sum(sum(slab.^2,2)./tausq);

        slab_temp = slice_sample(nslice, 0, log_dist_slab, slab, slice_width, false);
        slab = reshape(slab_temp(:,end),r,p);

        % sample spike by gibbs sampling
        for k = 1:(r*p),

            log_dist_spike = @(spike) sum(sum(Y.*(X*bsxfun(@times, slab, spike)') - log(1+exp(X*bsxfun(@times, slab, spike)'))));

            spike(k) = 0;
            mass0 = log_dist_spike(spike);

            spike(k) = 1;
            mass1 = log_dist_spike(spike);

            pr = 1./(1+exp(- (mass1 - mass0 + prior_odds(k))));
            spike(k) = single(rand(1,1,'single') < pr);
        end

        llik = sum(sum(Y.*(X*bsxfun(@times, slab, spike)') - log(1+exp(X*bsxfun(@times, slab, spike)'))));
        llik_samples(iter) = llik;

        if iter > burnin,
            spike_samples(:,:,iter-burnin) = spike(:,1:p);
            slab_samples(:,:,iter-burnin) = slab(:,1:p);
        end

        fprintf('iter = %05d, llik = %.4e\r',iter,llik);
    end


    effect = mean(spike_samples.*slab_samples,3);
    effect_var = mean((spike_samples.*slab_samples).^2,3) - effect.^2;
    effect_z = effect ./ max(1e-3,sqrt(effect_var));
    effect_p = 2*(1-normcdf(abs(effect_z)));

end
